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We study the spin dynamics of individual black holes in a binary system. In particular we focus on 
the polar precession of spins and the possibility of a complete flip of spins with respect to the orbital 
plane. We perform a full numerical simulation that displays these characteristics. We evolve equal 
mass binary spinning black holes for t = 20, 000M from an initial proper separation of d = 25 M 
down to merger after 48.5 orbits. We compute the gravitational radiation from this system and 
compare it to 3.5 post-Newtonian generated waveforms finding close agreement. We then further 
use 3.5 post-Newtonian evolutions to show the extension of this spin flip-flop phenomenon to unequal 
mass binaries. We also provide analytic expressions to approximate the maximum flip-flop angle 
and frequency in terms of the binary spins and mass ratio parameters at a given orbital radius. 
Finally we discuss the effect this spin flip-flop would have on accreting matter and other potential 
observational effects. 
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I. INTRODUCTION 

In this first decade since the breakthroughs HH20 that 
allowed to numerically solve General Relativity’s field 
equations for the evolution of black hole binaries (BHB), 
we have gained many new insights on these systems. The 
computation of theoretical gravitational waveforms from 
these BHB systems is very important for first detection 
and estimation of binary’s parameters 0HZ1- It is also 
of astrophysical interest to study the BHB orbital and 
spin dynamics in precessing systems and how accreting 
matter interacts with them. 

The spin of each individual black hole (BH) can no¬ 
tably affect its orbital motion as displayed, for instance, 
by the hangup mechanism [5], which delays or prompts 
the merger of the binary according to the sign of the spin- 
orbit coupling, thus ensuring the formation of a horizon 
(cosmic censorship hypothesis) at merger before the final 
hole settles to a Kerr BH nun] (no hair theorem). 

One of the most notable predictions of numerical 
relativity is that the remnant of the merger of two 
highly spinning BHs may receive a recoil of thousands 
of km/s |1H!T? due to asymmetrical emission of grav¬ 
itational radiation induced by the BH spins FI H5] , 
This developed intense astronomical searches for such 
fast moving BHs that may completely escape from their 
host galaxies or produce observable disturbances to the 
velocity field of stars in their cores [TB1J . 

New developments in numerical relativity allowed the 
study of relatively far separated binaries, up to distances 
of 100M in Ref. 17], and mass ratios of 100 : 1 in 
Refs. [18l fl9]. Very long term evolutions are now pos¬ 
sible mm as well as the study of near maximally spin¬ 
ning BHBs fUl - fUl 

In Ref. [20], we have found that the spin of BHs can 
completely reverse sign during its orbital stage. This 


flip-flop of spins is due to a spin-spin coupling effect and 
may have important observational consequences when 
the binary system is accreting gas from a galactic en¬ 
vironment 155]. In this paper we extend the analysis to 
unequal mass binaries using 3.5 post-Newtonian (PN) 
evolutions. Those proved a reliable description for this 
scenario when compared to the long term full numerical 
simulation of an equal mass binary [5DI . 


This paper is organized as follows. In Sec. [TT] we re¬ 
visit the equal mass binary scenario providing further 
detail and analysis of our full numerical simulation. We 
also provide a simple vector analysis to understand the 
mechanism behind the flip-flop and then a 2PN analytic 
study to describe its spin dynamics. In Sec. |III| we study 
the case of unequal mass binaries. We use 3.5PN evolu¬ 
tions of the equations of motion coupled with the 2PN 
spin dynamics to find the configurations that maximize 
the flip-flop angle and the likelihood of a flip-flop an¬ 
gle to occur given plausible astrophysical distribution of 
binary parameters. We also provide 2PN analytic expres¬ 
sions to approximate the flip-flop angle and frequency as 
measured in the orbital frame and asymptotic frame de¬ 
scribed in terms of projections with respect to the orbital 
angular momentum L and total angular momentum J. 
We end the paper with Sec. |IV| where we make some es¬ 
timates of the effect this flip-flop phenomena could have 
in realistic astrophysical scenarios and suggest that de¬ 
tailed simulations involving the magnetohydrodynamical 
description of accreting matter onto BHs are needed to 
find the characteristic electromagnetic signatures of the 
flip-flop. The appendix [A] provides details on the PN 
computation of the flip-flop frequencies and angles. 
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TABLE I. Initial data parameters and system details. The 
punctures are located at ri = (xi, 0, z) and fh = (* 2 , 0 , z), 
with momenta P = ±(0, P, 0), spins Si = (0, 0, Siz) 
and S 2 = {Six, 0, Siz), mass parameters m p , horizon 
(Christodoulou) masses m H , total ADM mass Madia, and 
dimensionless spins a = a/rriH = S/m 2 H . The horizon masses 
and spins are given after the gauge settles, and the errors in 
those quantities, denoted as Sm H and 5a, are determined by 
the drift in the quantity during the inspiral. Also provided 
are the simple proper distance d, eccentricity at the start of 
the inspiral a, and eccentricity e/ and the number of orbits 
N just before merger. 


xi/m 

Xi/m 

z/m 

P/m 

d/m 

10.73983 

-10.76016 

-0.01968 

0.05909 

25.37 

m\/m 

m p /m 

Si z /m z 

Si x /m 2 

Siz/m 2 

0.48543 

0.30697 

0.05 

0.19365 

-0.05 

MADM/m 

Jadm/w 2 

e. 

e/ 

N 

0.99472 

1.2704344 

0.0322 

0.0006 
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ml? /m 
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FIG. 1. Initial configuration of the simulated binary. 


II. EQUAL MASSES BINARIES 

We first study the case of equal mass binaries since 
they cleanly display the flip-flop effect on the spins of 
the holes and allow for a simple approximate analytic 
model of the process. In Sec. [Ill] we will study in detail 
its mass ratio dependence. 


A. Full Numerical Evolution 

In order to verify the realization of the spin flip-flops in 
the presence of gravitational radiation and the full non- 
linearities of General Relativity during the final stages of 
the inspiral, we consider a long-term full numerical sim¬ 
ulation with initial configuration as described in Table [I] 
and depicted in Fig. [I] 

We evolve the following BHB data sets using the 
LazEv [55] implementation of the moving puncture ap¬ 
proach Ell2| with the conformal function W = yfx = 


exp(— 2<j>) suggested by Ref. [27]. For the run presented 
here, we use centered, eighth-order finite differencing in 
space [25] and a fourth-order Runge Kutta time integra¬ 
tor (note that we do not upwind the advection terms). 

Our code uses the EinsteinToolkit j25J [30] / Cac¬ 
tus m / Carpet 32] infrastructure. The Carpet 
mesh refinement driver provides a “moving boxes” style 
of mesh refinement. In this approach, refined grids of 
fixed size are arranged about the coordinate centers of 
both holes. The Carpet code then moves these fine 
grids about the computational domain by following the 
trajectories of the two BHs. 

We use AHFinderDirect [33] to locate apparent 
horizons. We measure the magnitude of the horizon 
spin using the isolated horizon (IH) algorithm detailed 
in Ref. [34] and as implemented in Ref. [35j. Note that 
once we have the horizon spin, we can calculate the hori¬ 
zon mass via the Christodoulou formula 

% = \J toL + sy(4m? r ), (1) 

where Wi rr = \J A/{ 167t), A is the surface area of the 
horizon, and Sh is the spin angular momentum of the 
BH (in units of M 2 ). In the tables below, we use the 
variation in the measured horizon irreducible mass and 
spin during the simulation as a measure of the error in 
computing these quantities. We measure radiated energy, 
linear momentum, and angular momentum, in terms of 
the radiative Weyl Scalar -04, using the formulas provided 
in Refs. mm- However, rather than using the full ?/q, 
we decompose it into £ and m modes and solve for the 
radiated linear momentum, dropping terms with £ > 6. 
The formulas in Refs. [36], .37, are valid at r = oo. We 
extract the radiated energy-momentum at finite radius 
and extrapolate to r = oo, and find that the new pertur¬ 
bative extrapolation described in Ref. 38] provides the 
most accurate waveforms. While the difference of fitting 
both linear and quadratic extrapolations provides an in¬ 
dependent measure of the error. 

Figure [2] highlights the important effects of precession 
of the orbital plane during the whole evolution and the 
three precession cycles occurring over the nearly 50 or¬ 
bits. 

Through evolution we track the horizon mass and mag¬ 
nitude and components of the individual spins of the 
BHs. Figure [3] displays the levels at which these quanti¬ 
ties are conserved during the t = 20, 000M of full numer¬ 
ical evolution and this provides a measure of the accu¬ 
racy of the simulation, i.e., the variations in Amf/raf ~ 
4 x 10 -5 , AS? /S? « 3 x 10" 3 , A</m? « 1.2 x 1CT 5 , 
and ASf/Sf » 9 x 1CT 4 . 

The spin components in the simulation coordinate sys¬ 
tem are displayed in Fig. [4] They show the polar compo¬ 
nent of Si from alignment at t = 0 to almost misalign¬ 
ment at merger and the complementary behavior for S 2 . 
From the variations with time of the {x, y) components 
of the spins one can also read-off the precessional effects. 
While these are coordinate based periods, one can verify, 
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FIG. 2. Precession of the orbital plane as displayed by the 
distance vector d = x±(t) — X 2 (t). 
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FIG. 3. The conservation of the individual masses and spin 
magnitudes of the holes during the binary evolution. 


for instance from the oscillations in the amplitude of the 
(£ = 2, m = 1) mode of the gravitational waveform in a 
gauge invariant way that those periods correspond to the 
physical dynamics of spins (see Fig. [6] in this paper and 
also Ref. [59]). 

Another useful representation of the spins m is given 
in Fig. [ 5 ] where the flip-flopping spin’s (Si) trajectory can 
be described as a composition of the polar (flip-flop) and 
azimuthal (precession) dynamics which resembles that of 
an orange peeling. For the sake of future interpretations 
and applications we also provide the spin trajectories as 
seen in an (approximately) inertial frame where the direc¬ 
tion of the total angular momentum J is conserved and 
coincident (approximately) with the coordinate direction 
z. 

The gravitational waveforms from spinning, process¬ 
ing BHB systems are of great interest for detection 
and parameter estimation in the forthcoming observa¬ 
tion of gravitational waves by Advanced LIGO [3D] and 




FIG. 4. BH (1,2) spin components in the coordinate frame. 
The top panel is BH1 and the bottom panel is BH2. 




FIG. 5. Change in the spin direction (in blue) of the BHl 
(with the smaller spin magnitude) in the orbital frame, L 
(left), and the coordinate frame z (right). Plotted also the 
directions of L in red and J in green. 


other laser interferometric gravitational wave observato¬ 
ries. Our simulation produced one of the longest wave¬ 
forms studied to date, starting from an initial configura¬ 
tion with proper separation, d = 25 M. This is well in 
the post-Newtonian regime (see also the recent simula¬ 
tion |2lj for nonspinning, unequal mass binaries starting 
at an initial separation d = 27 M). In Fig. [fl] the leading 
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FIG. 6. The full numerical (£ = 2, m = 2) waveform (above) and (£ = 2, m = 1) waveform (below) as extracted by observers 
at r = 175 M. 




modes (£ = 2, m = 2) and (2, 1) are displayed i2tjj and 
the first shows the well known chirp behavior with es¬ 
sentially monotonic increase of amplitude and frequency 
until merger. The second waveform, corresponding to the 
(£ = 2, m = 1) mode, displays a modulation in ampli- 


Of particular interest is to determine the precision of 
this waveform. Due to the high computational demand of 
this simulation, 2.5 million service units on 25 to 30 nodes 
of our local cluster “Blue Sky” with dual Intel Xeon E5- 
2680 processors nearing 100M of evolution per day, we 
have not produced detailed convergence studies involving 
several runs at different resolutions. Instead we compare 
our full numerical waveforms with those produced with 
3.5PN evolutions and found excellent agreement, partic¬ 
ularly at early times, when 3.5PN is expected to be very 
accurate, given the large initial separation of the holes. 

Figures [7] and [8] display the differences of the numer¬ 
ical and PN waveform for the (£ = 2, m = 2) mode 
phase and amplitude, respectively. Since the full numer¬ 


tude that corresponds to the period of precession of the 
orbital plane (and in this equal mass binary case this is 
coincident with the precession of the spins). Thus this 
waveform provides an independent and gauge invariant 
measure of the precession period. 


ical waveforms are extracted at a finite observer loca¬ 
tion (in this case at r = 175M) we extrapolate them to 
null infinity with a recent 0(l/r 2 ) accurate perturbative 
method Ei E] to obtain further agreement, showing 
reasonable errors until near the final merger. 

Since this BHB simulation is relatively long-term it is 
well suited to study the evolution of the eccentricity as 
a function of the binary’s separation. Using the method 
of Ref. [52] to determine e^, we measure the eccentric¬ 
ity per orbit (we have not included the first two orbits, 
contaminated by gauge settling, and last few orbits due 
to strong radiative effects). As the binary’s separation 
shrinks due to gravitational radiation the eccentricity re¬ 
duces as shown in Fig. [9j We performed a fit to the mea- 
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FIG. 7. The phase difference between the full numerically 
generated (t. = 2, m = 2) waveform (NR) and the 3.5PN 
waveform. Differences are notably reduced when the NR- 
waveform, extracted at r = 175A/, is extrapolated to an in¬ 
finite observer location by the perturbative formulae (PE) of 
1st order, accurate to (1/r) terms and second order, accurate 
to (1/r 2 ). 


CM 


E 

< 



FIG. 8. The amplitude difference between the full numer¬ 
ically generated (£ = 2, m = 2) waveform (NR) and the 
3.5PN waveform. Differences are notably reduced when the 
NR-waveform, extracted at r = 175 M, is extrapolated to an 
infinite observer location by the perturbative formulae (PE) 
of 1st order, (1/r), and second order, (1/r 2 ). 


sured eccentricity versus separation of the form ar b ob¬ 
taining an exponent b = 1.73486±0.1495 to compare with 
the lowest post-Newtonian prediction of Ref. [43] for ra¬ 
diation of eccentricity, e ~ r 19 / 12 . With 19/12 = 1.5833, 
the post-Newtonian prediction is within the statistical 
error bars of the fit. 

After merger the binary forms a single final remnant 
BH with characteristics given in Table [H] Notably the 
recoil velocity reaches 1500 km/s in agreement with the 
empirical predictions nmm. 

Figure [To] displays the settling of the final mass and 
spin measures as the final BH radiates away its last dis- 



FIG. 9. Evolution of the eccentricity versus coordinate sep¬ 
aration of the BHB. Dots represent measurements of the ec¬ 
centricity, en, and the continuous curve a fit to its decay with 
e ~ r .i-73486±o. 1495 j n com p ar i son] the theoretical prediction 
is e ~ r 1 ' 5833 . 


TABLE II. Remnant properties and recoil velocity. The fi¬ 
nal mass and spin are measured from the horizon, and the 
recoil velocity is calculated from the gravitational waveforms. 
The error in the mass and spin is determined by the drift in 
those quantities after the remnant settles down. The error in 
the recoil velocity is estimated by the difference between first 
and second order polynomial extrapolation to infinity of the 
waveforms. 


M-rem / TYl 

| (%rem \ 

Vrecoii [km/s] 

0.94904 ±0.00000 

0.70377 ± 0.00002 

1508.49 ± 16.08 

OL X 

rem 

a y 

rem 

Oi Z 

'- Aj rem 

0.10815 ±0.00003 

-0.01986 ±0.00000 

0.69513 ± 0.00002 


tortions and becomes a quiet Kerr BH [5]. 


B. Vector Analysis 

In order to visualize the origin of this flip-flop effect 
we can use a simple vector addition model to represent 
the total angular momentum of the system, J = L ± S , 
as the sum of the orbital angular momentum L and the 
total spin S = Si + S 2 . This total spin being the sum of 
the individual spin of the holes Si and S 2 as represented 
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FIG. 10. The mass and spin of the final merged BH. 
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FIG. 11. Flip-flop of Si and S 2 as they precess with L. 


in Fig. [TT] 

In this picture, the six degrees of freedom of the spins 
Si and S 2 are shown as vectors in 3-space. These de¬ 
grees of freedom are constrained by the conservation of 
the magnitude of the individual spins Si, and S 2 , the 
magnitude of the total spin S, and its projection along 


the orbital angular momentum, S • L. The conservation 
of the magnitudes of the individual spins and the mag¬ 
nitude of their sum S, can be seen in the 2PN evolution 
of the spins in Eqs. |7| of Sec. IIC The conservation of 
So • L was shown in Ref. [25] using the orbit averaged 
version of Eqs. ([7]). The conservation of these quantities 
in full numerical simulations of BHBs has been observed 
to be approximately true in Refs. [HI] and m- Note 
that in what follows to determine the flip-flop angle, we 
do not use the PN (or full numerical) dynamics beyond 
the referred conservation of the individual and total spin 
magnitudes. 

It follows hence that the following quantities are con¬ 
served (see Fig. 12 1 


S ■ S = S 2 = Sf + S% + 2 S 1 S 2 cos (3 = constant, 

S • Si = SSi cos 7 = Si + S 2 S 1 cos f3 = constant, 

S • S 2 = SS 2 cos (/3 - 7 ) = S 2 2 + S 2 Si cos /3 

= constant. ( 2 ) 


FIG. 12. The configuration of the spins in the orbiting frame. 
Spin configurations Si and S 2 relative to the orbital angular 
momentum L. Here S = Si + S 2 . 


Thus both spins, Si and S 2 , oscillate around S which in 
turn precesses around L. 

The direction of oscillation of the spins can be seen 
from the leading precession equation for dL/dt = 
(7/2r 3 )(J x L ) indicating that L oscillates counterclock¬ 
wise. Likewise the leading precession for dS\/dt = 
(l/r 3 )((7/2) J x Si — 3S x Si) indicates that Si oscillates 
clockwise around S (see Eqs. © for the precise expres¬ 
sion). Thus Fig. [TTJ is represented as seen by an observer 
located up (north) with respect to J direction. 

We can define the flip-flop angle with respect to the po¬ 
lar coordinates as measured from L (analogous formulae 
can be defined for J) 

A 8{ f = d™ ax - 9™ in , (5) 

where i = 1,2 labels either BH spin. 

We have three cases to consider depending on 63 , 
where cos ( 63 ) = S^/S = (S ■ L)/S, 

( 9 s + 7 if 0 < 0 S < 7 , 

A 9{ f = < 27 if 7 < 9 s < 7t — 7, (6) 

I 7 T — 63 T 7 if 7 T — 7 < < 7 T , 


This leads to the conservation of (3 and 7 during the 
evolution of the binary. 

In particular we find that <ST oscillates around S be¬ 
tween angles 7 and —7 (when it is both coplanar to S 
and L ), where 

S 1 + S 2 cosP _ S 2 + S 2 - S 2 ^ 
C0S 7 ~~ yjS 2 + S 2 + 2S 1 S 2 cos /3 ~~ 2SSi 

BH2 also oscillates at the same flip-flop frequency tiff, 
but with a smaller angle (since we consider S 2 > Si) 
given by ±(/3 — 7 ) where 


cos (/3 — 7 ) 


S 2 + S 2 - S 2 
2SS 2 


(4) 


for 0 < 7 < 7 t/2 , and similar for BH2, replacing 7 by 
/3 — 7 . When 7 r /2 < 7 < 7 r, we replace 7 f> tt - 7 in the 
equation above. 

Hence, for instance, in order to maximize the flip-flop 
angle, we can set 7 = 7r/2, which leads to the vanishing 
of cos 7 on the left hand side of the above equation ©, 
leading to the condition Si + S 2 cos (3 = 0, for q = 1, 
used in our full numerical configuration, as detailed in 
Sec. UTAl 

This oscillation of the spins represent a genuine spin- 
flip in the sense that it is the same object that completely 
changes its spin orientation. This is different from the 
simple case where the final remnant spin has flipped di¬ 
rection compared to the spin of one of the individual or¬ 
biting BHs [35]. It is also different from a classic analog 
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of the hyperfine transition of an electron in the Hydrogen 
atom flipping its spin and leading to the famous 21cm 
radio emission. Here the flip also occurs in a conservative 
set up, with no emission of gravitational waves, and it is 
compensated by an opposite flop from the other BH. 


C. Post-Newtonian Spin Dynamics 


The vector analysis is helpful to visualize the angles of 
flip-flop but in order to determine their rate of oscillation 
one has to resort to a dynamical computation. In order to 
provide simple, approximate, expressions we will perform 
a 2PN (conservative) study. 

The precession equations for the spins Si and S2 with a 
mass ratio q = mi/m 2 to leading spin-orbit and spin-spin 
couplings in the (2PN) post-Newtonian expansion [48] 
take the form 


dSi 

dt 

d§2 

dt 


1 

/3 

1 

y3 



L — S -2 + 


3(S 0 • n) „ 

— - n 

1 + q 


xSi, 



L — §1 + 


3 q{S 0 ■ n) A 

—7 - n 

1 + Q 


x S 2 ,(7) 


where n = (rq — r 2 )/|fq — r 2 ] and 


for Si and similar equations for S 2 by exchanging labels 1 
and 2 from Eqs. Q . In the above expression, v\ = r£l or t, 
is the tangential velocity of the binary (in a quasicircular 


orbit, see Eq. (131 for Q orb ), SA = <S) • L, S ifl = S* • n, 


S iX = Si ■ A and S L = S ■ L. 

We can then obtain equations of the form d 2 (Si ■ 
L)/dt 2 = Q 2 j: Si ■ L + ■ ■ ■ for i = 1, 2 and analogously 
for the perpendicular component of Si giving From 
where we can read-off the orbit averaged polar and az¬ 
imuthal oscillations frequencies of the spin Si (see also 
Ref. [IH for the equations projected along J) 


Vtff = 3 


S 


1 - 


2 S-L 

M 3/2y.l/2 


79 9 ^ - 

n P = ^3+ S ' L ) + 


( 11 ) 

( 12 ) 


that we identify with the flip-flop and precession frequen¬ 
cies respectively. For the sake of completeness and fur¬ 
ther reference, the orbital frequency of equal mass bina¬ 
ries in quasicircular orbits is given by Ref. [51] 


MQ, or b 


/m \ 3/2 11 /m \ 5/2 

It ) -1(7) + 


(13) 


So — ( 1 H— ) Si + (1 + q) S2 . 


( 8 ) 


For more details see, for instance, the reviews in Refs. [491 
150] , The averaging of Eqs. ([7]) over the orbital period 
gives 


(A) = ^ 


<&> = 


1 


3^ 

2 t 

3 St 

2 l 

3 St 

2 t 

3 S t 


L + \S2 

J-3S 


L+ 2 Sl 


= — I [ — — - -jf- } J — 3S 


x Si 
x Si, 
x S 2 
x S 2 , 


(9) 


where £ = \L\ = qM z l 2 r x l 2 f{\ + q ) 2 and we have ignored 
cubic terms of spins in the most right hand side of each 
equation. 

By decomposing the spins along L and perpen dicu lar 
to it with unit vectors A and h, as shown in Fig. 
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the fashion of Sec. IV.A of Ref. [44], we obtain the spin 
evolution equations for q = 1 as 


Sin = S lX 


S lX = -S lft 


Sli 




+ 


Si L = 


!*_".( i + 
r 2r 3 \ ' 

"-(1 + 
r 2 r 3 1 

I 9 S 1 iS 2 n _ 

' 2 r 3 r 3 

9 S lft S 2 ^ 9S 1 ^S 2ft 3S^S, 


S 2 l 


9SilS 2 \ 

2 r 3 


S 2 l 


3S 1 fSf, 




2 r 3 


2 r 3 


1+^1, ( 10 ) 


In order to verify the accuracy of the above expressions 
for the flip-flop frequency and angle we performed a series 
of 3.5PN numerical integrations representing the merger 
of BHBs from r = 100M down to r = 5 M and measured 
the flip-flop angle and the frequency of oscillation. The 
results are presented in Figs. [13| and [14| as a color map 
with the direct integration of the PN equations of motion 
coupled to the spin evolutions and we superposed a few 
curve levels based on the above analytic estimates. We 
observe good agreement among them, in particular for 
the maximum flip angle described by cos/3 = —Si/S 2 - 
Fig. [l4| displays as a color map the frequency at the start 
of the evolution, r = 100M where the term (3 S/r 3 ) in 
Eq. ( [IT] ) dominates. In this figure, the thin lines are the 
contours from the measured flip-flop frequency, and the 
bold lines are the analytic expression given by Eq. © 
for quasicircular orbits at the corresponding r. 


D. Statistics 

Let us study first the simplest case of conservative pre¬ 
cession for equal mass binaries. In Fig.|15[ we display the 
probability of a flip-flop angle A 9ff>x assuming equal 
masses with random spin orientations and magnitudes of 
BH1 and 2. For the sake of simplicity, we used here spin 
evolution equations 0 with no radiation reaction taken 
into account. The line labeled ‘2D’ considers random an¬ 
gular distribution of one spin and random distribution of 
the ratio of spin magnitudes. It could represent the sce¬ 
nario where accretion processes succeeded in aligning the 
spin of one of the BHs, while the spin of the other has still 























by 


Maximum flip angle 

180 

160 

140 

120 

100 
CO. 

80 
60 
40 
20 
0 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

S/S 2 


FIG. 13. q = 1 binaries, inspiraling from r = 100M to r = 
5 M. Initial configurations having (pi=0 = (f> 0? = 0, and 
varying 9° = /3 and S 1 /S 2 ■ The central blue curve follows 
cos f3 = —S 2 /S 1 and agrees with the maximum (yellow) of 
the flip-flop angle of 62 - The other two blue lines represent 
the 90-degrees flips level. 
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FIG. 14. Flip-flop frequency oscillation for q = 1 binaries, 
inspiraling from r = lOOAf to r = 5 AT Starting from initial 
(j>i = 0 = (j>%, 9 1 =0, for different choices of 9® = /3 and 
Si/S 2 - Different surface levels are represented in units of 
10 - 7 for the initial frequency at r = lOOA-f. Thin lines are the 
measured flip-flop frequency during evolution and thick lines 
their analytic approximation for quasicircular orbits, Eq. (111. 


a random orientation due to a much larger time scale for 
alignment. If we allow instead to vary all six parameters 
of the BH spins (magnitude and both spin directions) 
with q = 1, we observe the distribution labeled ‘ 6 D’. 

When either BH has its spin aligned or counteraligned 
with the orbital angular momentum and the other spin 
is chosen randomly we have a random angle /3 between 
the directions of the BH1 and 2 spins. In this case one 
can obtain from Eq. © that the probability of an angle 
7 given a random distribution of cos /3 is approximated 


P(cos 7 ) = 1 + (2cos7 - 1) , (14) 

which upon assuming random distribution of spin mag¬ 
nitudes 0 < (S 1 /S 2 ) < 1, leads to the probability for a 
flip-flop angle larger than x 

P(Ad ff > x) = i cos (cos (0 + l) . (15) 

This probability distribution, represented in Fig.|l5|by 
the red continuous curve shows, for instance, that (for 
equal mass binaries) there is nearly a 60% probability of 
flip-flops of more than 90°. 

In the more general case, allowing for radiation re¬ 
action to drive the BHBs towards merger and random 
variation of all variables, we find the same qualitative 
behavior as above plus a detailed A0-dependence. From 
the expressions for the flip-flop angle in Eq. |6]), we ob¬ 
serve they depend entirely on the two angles 7 and 9s, 
both being constants of motion (for equal mass binaries). 
Thus once we express them in terms of initial data we can 
predict the flip-flop angle. This can be achieved in terms 
of four of the six spin components: 9\, 62 , A (f> = <^2 — <t>i, 
and Si/S 2 , from Eq. ([ 3 ]) for 7 

*^ 1/^2 T cos / 

cos 7 =— . , (16) 

a/1 + Sf/S% + 2 S 1 /S 2 cos/3 

and 

S 1 /S 2 cos 0 \ + cos 0 2 

cos 8 s = , , (17) 

^l + Sf/S'l + 2 Si/S 2 Cosp 


where 


cos /3 = cos 8 \ cos 82 + sin 9\ sin d 2 cos Ac\>. 


(18) 


The largest flip-flop angles occurs when the two BH 
spins started at (or pass through) some initial configu¬ 


ration with A </> = 7t/ 2 as is displayed in the Fig. [16 
The curves in this figure represent the probability that 
an equal-mass system with initial A cj> will have a flip-flop 
angle > x given random spin orientations and spin mag¬ 
nitudes of the primary and secondary BHs (here labeled 
simply as 1 or 2). Each curve corresponds to a set of 
203,401 PN evolutions (|ai| x |o :2 1 x cos(0°) x cos^) 
= 11 x 11 x 41 x 41) starting from a separation of 
r = 100M and evolved with 3.5PN radiation reaction 
terms down to a separation of 5 M. The color of the 
curve is determined by the initial value of A (f>. 


III. UNEQUAL MASSES BINARIES 

The flip-flop of BH spins in a binary system is an in¬ 
teresting effect since it singles out a conservative dynam¬ 
ical behavior that may produce important astrophysical 
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FIG. 15. The probability for the q = 1 case of having a flip 
greater than a certain angle. Distribution ‘2D’ lets either 9i 
or #2 be aligned or counteraligned with L with a randomly 
varied cos of the other variable and spin magnitudes ratio. 
There are 209,450 samples that match the analytic distribu¬ 
tion. Distribution ‘6D’ shows the probability of flip over a 
generic range of spins (6 parameters) for q = 1, totaling 16 
million evolutions. 
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Flip Angle (degrees) 

FIG. 16. The probability of a flip-flop angle A 9;f > x assum¬ 
ing equal masses, (q = 1), and a random initial spin orien¬ 
tations and magnitudes of both, the primary and secondary 
spins. Evolutions from initial separation r = 100M down to 
r = 571 1 are performed using 3.5PN equations of motion. The 
dependence with the azimuthal orientation A<j> = d(f> is dis¬ 
played. We observe the largest flip angles probabilities for 
configurations with A <j> = 7r/2 and smallest for A(p = 0 and 

7 T . 


phenomena. In order to further evaluate the relevance to 
astrophysics we have to study BHBs with unequal masses 
and the likelihood of a given flip-angle to occur as well 
as the timescale (or flip-flop frequency) involved. 

In the equal mass case, we have been able to discuss 
the geometry of the flip-flop angles in Sec. |IIB| based on 
(to a good approximation) the conservation of the spin 
magnitudes, Si, S 2 , and S as well as the component of 


the total spin along the orbital angular momentum S^. 
These four conserved quantities allowed us to express the 
flip-flop angle in terms of the conserved quantities 7 and 
0 S , i.e., Eq. <©• The conservation of 7 and 0s implies 
that we are able to determine the flip-flop angle in terms 
of the initial configuration of the binary, and hence by 
specifying S 1 /S 2 , 0 1 , 0 2 , and A <j>. Note that from the six 
variables that are needed to specify the components of 
the two vectors Si and S 2 only four are needed here due 
to the fact we compute a dimensionless flip-flop angle, 
hence the dependence on the ratio S\/S 2 and the depen¬ 
dence on A<t> = (f> 2 — (pi appears due to the fact that 
for large separations a rotation of the angle <f> is already 
provided by the precession of S around L, spanning all 
orientations relative to the orbital linear momentum, p). 

In retrospective, while the flip-flop angle for equal 
masses can be determined geometrically, the frequency at 
which spins flip-flop is a dynamical quantity that requires 
information about the evolution of the spins as well as the 
orbital evolution of the BHBs. At the lowest (2PN) post- 
Newtonian approximation we can evaluate this frequency 
using the spin-evolution equations (|t| assuming that the 
binary is separated far enough that the radial decay is 
negligible during a flip-flop cycle (note that the gravita¬ 
tional radiation frequency scales as ~ r -4 compared to 
the flip-flop one scaling as ~ r -3 ). The flip-flop frequency 
thus depends on the orbital radius r, the total spin and 
its projection along L as given by expression (111]). 

Here, in the unequal mass binary case we can assume 
to a certain degree of accuracy the conservation of both 
spin magnitudes, Si and S 2 . S £ is not conserved, but the 
projection of the vector So (as given in Eq. (|sb) along L , 
denoted as S Q ^ is approximately conserved H5| . How¬ 
ever, the total spin magnitude, S, is not conserved since 
the angle between the two individual spins (3 is no longer 
conserved. We will use a combination of analytic and 
3.5PN numerical integrations to provide asymptotic ex¬ 
pressions for the flip-flop angle, flip-flop frequency, and 
estimates of the likelihood for them to arise in astro- 
physical scenarios. Contemporary studies [52| could also 
describe the spin flips in terms of the two extreme values 
of the total spin. 

Note that in generic binaries spin resonances of the 
kind studied in Refs. [53 M bring the azimuthal angle 
differences towards A<j> = 0 and 7 r when the spin polar 
orientations are significantly different, and to A cj> = 7t/2 

when they are similar [55H5Z1; thus allowing a further 

statistical reduction of the parameter space to explore. 


A. Post-Newtonian Analysis 

In the unequal mass binaries case the flip-flop angle will 
be a function of not only the intrinsic spins, a, = S t /rrif, 
but also the radius r of the quasicircular orbit as well as 
the mass ratio q. 

By studying the spin evolution equations ([T]) projected 
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along L (or along J) we can obtain the maximum flip-flop 
angle for a given set of binary parameters, oq, 012 , q, and 

r, 


1 — cos(A${^) 


2a! 


(1 -q¥ 


M ^ | 4aia!g fM 
(1 - q) s 


3/2 


(19) 

This maximum flip-flop angle of the smaller BH with 
respect to L is achieved when the smaller hole starts 
(or pass during its oscillation) through an anti-alignment 
with respect to L, i.e., 9® L = n and the larger BH spin 
S 2 forms initially an angle with L given by 


cos 


O' 2 


3qaiCt2 f M 


( 20 ) 


On the other hand, if we seek to maximize the flip-flop 
angle of the larger BH, we find 

( 21 ) 

This maximum flip of the larger hole is achieved when the 
spin of this hole goes through alignment with the orbital 
angular momentum, i.e., = 0, then the smaller BH 

spin forms an angle ^ with S 2 given by 


cos 9° 


1 L 


qa 1 
(9-1) 


M (-) ■ <*» 


In Fig. [TT] we display the results of 3.5PN evolutions 
from r = 100M for 10 choices for q , 41 choices for cos 6 ®, 
41 choices for cos 02, 11 choices for |cki |, 11 choices for 
| Oi 2 1, and 37 choices for A<j> by setting </>° = 0 and choos¬ 
ing different (fy. The systems are then evolved with ra¬ 
diation reaction down to 5 M. The color-bar in the fig¬ 
ures denote the initial A<f>. The highest values happen at 
A (j> = 7t/ 2, but the results are very insensitive on this de¬ 
pendence for almost every mass ratio below q ss 0.7. This 
same property is observed in both the L- and J-frames. 

We verified these analytic dependencies by performing 
numerical integrations of the spin evolution equations 0 
for different radii, r/M = 30, 50, 100, 150, 200, 250, and 
mass ratios q = 1/2, 2/5, 1/3, 1/4, 1/5 and display a 
summary in Figs. [18] and [l9j We compare in those 
plots (not fit) the analytic form of the leading flip an¬ 
gle (on the right-hand side of Eqs. (191 and ©) ver¬ 
sus the measured values as obtained from assuming the 
forms l-cos(6 l {{) = /r + B Li s /r 1 ’ 5 + Cl,s A" 2 , and 
= DL,s/r + EL,s/r 15 + FL,s/r 2 - In order to make 
explicit the 1 o 2, g <-»• 1/q symmetry of the A and D 
coefficients, we display the case of the larger BH flip-flop 
on the left of Figs. 18 and[l9|in terms of 1/q > 1 variable. 


The maximum angle of flip-flop with respect to J- 
frames is different than with respect to the L-frame due 
to the fact that, even during conservative evolution, the 
angle of precession of L around J is not preserved. In 




Flip Angle (degrees) 
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FIG. 17. The dependence of the probabilities of a flip-flop 
larger than a given value as a function of the initial A(f> an¬ 
gle among spins in the L-frame. Evolutions from r = 100 M 
down to r = 5 M are performed with 3.5PN including ra¬ 
diation reaction. From left to right and top to bottom for 
q = 0.9, 0.8, 0.7, 0.6. We observe an insensitivity on A <j> for 
mass ratios of the binary below q = 0.7. Similar results are 
observed in the J-frame. 




FIG. 18. Measured maximum flip-flop angle from PN evo¬ 
lutions versus analytic prediction (no f ittin g her e). Here we 
display the leading coefficient of Eqs. ( |19| ) and ( |21[ ) for the 
larger BH flip-flop angle (left) and smaller BH angle (right). 


fact, because for unequal masses, unlike the equal mass 
case, the projection of the total spin S along Z, S^, is 
not preserved, 


J • L = 


(e+s L ) 


(23) 


even if on average the magnitudes of the total angular 
momentum, j, and the orbital angular momentum, £, are 
preserved. 

This leads to the possibility of having larger fluctua¬ 
tions in the angles as “seen” in the J-frame, given by 


1 — cos(A 0[j) 


2a 2 2 (M\ 

q( l - q ) 2 V r ) 


and 


1 - cos(A 0 f 2 f j) 


2 afq 3 / M\ 

W 
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. towtfP,- 1 ) /«V' ! 

(9-1)3 V r / 

for the maximum smaller and larger BH oscillations of 
their spins respectively. The extreme values happen at a 
initial coordinate A cf> = 7r/2 and seem to correspond to 
the S ma x (and S m i n ) configurations in Fig. 2 of Ref. j58] . 
Note the q —> \/q symmetry of all the above expres¬ 
sions corresponding to a 1 —> 2 exchange in the labels of 
the holes. Also at each moment of time, the conservation 
of Sq determines the angular location of the other BH 
given the spin of the BH with larger flip-flop angle. 

The initial angle of the small hole (and between the 
spins) that maximizes the flip-flop as seen in the J-frame 
is approximated by 


COS0°j 


1 (1 + q) qa.\ 

'2 (1-g) 

1 ( 9 2 + 49 + l)aia 2 ( M \ 


2 (I-9) 2 

for the larger BH spin initially pointing to 


(26) 


cos e\j 


1 1 2 2 

1 - 2 qa ' 


qa i 2 a 2 


3/2 


(27) 


While if we seek to maximize the flip-flop angle of the 
smaller hole as seen in the J-frame, we find 


cos Oij 


2 9 2 


a 2 2 a\ 


3/2 


~) ( 28 ) 


In this case the spin of the larger hole points to 

Cm 


cos 0° 


1 (1 + 9 ) a 2 

27 ~ 2 9(1-9) 

1 (q 2 + 4q + l)a 2 ai f 

+ 2 (l^F 


- ■ (29) 


r ) 


The detailed derivation is shown in Appendix [A] 

The frequency of the flip-flop oscillation in all these 
cases is given by (see Appendix p\j) 


MCi{{ 


3 11 — 9 I (M \ 5/2 


2 (1 + 9 ) V r ) 


+3 


St-S ; 

M 2 


L ( M ' 

2 sign(l — 9 ). (30) 




FIG. 19. Measured flip-flop frequency versus analytic predic¬ 
tion (no fitting here). Here we display the leading coefficient 
of Eq. (301 for the larger BH flip-flop angle (left) and smaller 
BH angle (right). 


2a£(l + 49 ) + 90^(3 + 2 9 + 5g 2 ) 
2(1+ 9) 3 



The leading terms of the precession or azimuthal fre¬ 
quency are fl5j 


A 1 ' ■— 1 

+ 


7 9 /M \ 5/2 

2 (1 + q ) 2 \ r ) 

'7 {a[q 2 + a%) _ 3 {a[q + a%) 
.2 (I + 9) 2 4 (I + 9) 



If we now request that the flip-flop angle be exactly 7 r in 
an unequal mass configuration, this might only happen 
for quasicircular orbits inside a critical separation re¬ 
in order to estimate this critical separation we may use 
the analytic expression for the maximum flip-flop angle 
of the smaller BH given in Eq. (19) and request that 
A 0 {[ = 7 r for total flip-flop in the L-frame. This leads to 
the condition 


1 = u 2 + 2 pit 3 , 


(33) 


where 


u = u L = 


a 2 M L qai 

- 7 \ — . P = V = -■ (34) 

(1 - 9) V r c a 2 


The exact solution to Eq. (33) is given by 
1 

6 p 

where 


” = i (/ ,/3 + 


(35) 


The origin of the additional term for unequal masses 
scaling with ~ j- -5 / 2 is due to the non-conservation of 
the angle /3 between the two spins (as opposed to its 
conservation in the 9 = 1 case). These oscillations in /3 
are due to the differential precessional angular velocity of 
Si and S 2 for 9 ^ 1 and hence provides the (precessional) 
scaling r -5 / 2 . 

For the sake of completeness and comparison of the 
time-scales involved, we give here the leading terms of 
the orbital frequency in a PN expansion [5§ : 


Mri or b 



3 + 59 + 3q 2 (M \ 5/2 
2(1 + q) 2 W 


/ = (Vl08x/27p 2 - 1 p+ 54p 2 - l) . (36) 


It is useful though to estimate the best scenario for 
large rc, which occurs in the small to intermediate p- 
regime where u ~ 1. Thus we obtain 


Ma\ 


(37) 


with a best scenario for highly spinning BHs we get that 
9 > 3/4 is needed for rc > 20 M. Setting up a flip of 
90-degrees moves the critical radii outwards by approx¬ 
imately a factor two. This indicates that large flip-flop 
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angles, as measured in the orbital frame, are mostly ex¬ 
pected for comparable mass binaries. 

Similarly, to achieve the total flip-flop angle as seen in 
the J-frame, the above analysis applies with 


u = Uj = 


«2 


(1-9) 



j ai (3 -q)q 3/2 
P = P = 


2«2 


( 38 ) 

This leads to similar bounds as above for comparable 
masses, but gets a notable increase for small mass ratios 
due to large precessional effects (see the bottom panel of 
Fig. m. 


B. Statistics 


In order to evaluate the probabilities of a flip-flop an¬ 
gle to occur, we have performed numerical integrations of 
the 3.5PN equations of motion [251 EQi coupled with the 
spin evolution ([T]) for randomly varied spin magnitudes 
and orientations (thus representing unbiased initial con¬ 
ditions). This may correspond to ‘dry’ mergers of BHBs 
or to the case of ‘chaotic accretion ’ m■ Other scenarios 
involving spin alignment due to coherent accretion are 
possible [62] as well as the effect of resonances [53], that 
lead to two ‘attractor’ configurations, either A</ = 0 or 
A 4 > = 7T. We will study those two cases below. 

We have thus performed a set of ctq x cos(0[ ) ) x a 2 x 
cos($ 2 ) = 20 x 90 x 20 x 90 = 3, 240, 000 simulations per 
q (we have masked the points that have one of the spins 
0 or that have anti-aligned or aligned spins bringing the 
number for the statistics down to 2, 859,120). 

Figure 20 displays the results of the probabilities (as 
inferred from the occurring frequency of each value) of a 
flip-flop angle larger than a given (by the abscissa) value 
to occur for each mass ratio q in the range 0.1 < q < 0.9 
(no resonances exists for the q = 1 case). We observe 
that for some values of q > 0.7 the probabilities for large 
angles can be larger than for q « 1 indicating that the 
leading term in Eq. © is preponderant as the binary in¬ 
spirals down to smaller separations. This is a phenomena 
that we discussed in the context of maximal flip-flop at 
the end of Sec. Ill A We thus conclude that large flip-flop 
angles are frequent in comparable masses binaries. 

To further illustrate the probabilities of a given flip- 
flop angle to occur, we consider a BHB system bearing 
a mass ratio q = 1/2 at a relatively short distance, r = 
20 M. Figure [2l] displays the results of a uniform initial 
distribution of spin magnitudes and angles. 

We observe that the distribution of a given flip-flop an¬ 
gle has a peak at about 41° corresponding to the smaller 
hole flip-flop and at around 17° for the flip-flop angle 
of the larger hole. Flip-flop angles larger than 63° are 
very unlikely, in agreement with the results displayed in 

Fig.m 

We mention that for unequal masses there is a differ¬ 
ence in counting flip-flop angles in the J-frame (total an¬ 
gular momentum) and the J-frame (orbital angular mo- 
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FIG. 20. The distributions of the flip angles larger than a 
given (abscissa) value for a binary evolved with 3.5PN from a 
separation r = 100M to r = 5 M for each 0.1 < q < 0.9 in the 
J-frame. The upper panel assumes an initial <ji l = 4> 2 while 
the lower panel assumes an initial </>i = —cj >2 as the resonant 
attractors at larger separations. 


mentum). In Fig. 22 we display the differences for a few 
selected mass ratios. We observe in general that larger 
angles are seen in the J-frame. Particularly for q = 0.1, 
due to the effects of transitional precession [53), that may 


reverse the orientation of J. But we also observe larger 
angles for q > 1/4, configurations that exclude total tran¬ 
sitional precession [44lf47l ). 

Both, the J- and J-frames are useful references to ac¬ 
count for different astrophysical accretion scenarios. 


IV. CONCLUSIONS AND DISCUSSION 


The unequal mass component of the flip-flop frequency, 
as given by Eq. (30), adds a leading term with a depen¬ 
dence on the orbital separation as ~ r -5 / 2 . This term, 
also proportional to (1 — q), eventually becomes the domi¬ 
nant time-scale at large enough separations. We can com¬ 
pare this time-scale with that for the realignment of BH 
spins via the Bardeen-Peterson torque effect produced by 
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FIG. 21. The distributions of the flip angles for a binary at 
a constant separation r = 20 M and mass ratio q = 1/2. Spin 
magnitudes and angles have been chosen uniformly and the 
max(small, large) flip angles is evaluated. The peak at around 
41° corresponds to the smaller hole spin flip. The larger hole 
gives an additional peak (at around 17°) towards low angles 
and bulks up the middle region between peaks. 


accreting matter. This alignment mechanism competes 
with the spin flip-flop one. 

The leading flip-flop period is now given by 


T ff 


2,000 yr 


(1 + g) 
(1-9) 



which is much shorter than the gravitational radiation 
scale [6T] 


T olv ~ 1.22 x 10V (jFm)' 


( M 


v 10 8 Mq 


(40) 


and hence can take over the Bardeen-Peterson alignment 
mechanism 164] at an order of magnitude further sepa¬ 
rations (i.e., 10, 000M) than the gravitational radiation 
decay, used as reference time-scale in Ref. [65]. Thus, the 
flip-flop of spins shortens the time available for accretion 
to completely align the spins of the BH with the orbital 
angular momentum. However, according to Eq. (19) the 


angle of flip-flops gets notably reduced at large distances, 
hence their statistical relevance can be focused to the lat¬ 
est stages of merger, i.e., at closer separations or compa¬ 
rable mass binaries. 

While in this paper we point out the relevance of study¬ 
ing accretion and torque exchange of matter with BHs 
in a binary system rather than on single BHs, a defini¬ 
tive account of these competing mechanisms should be 
provided by full BHB simulations in a matter filled en¬ 
vironment that incorporates all the needed elements of 
a magnetohydrodynamical description of the accretion, 
like in Ref. [66] with spinning BHB backgrounds. 

Since this flip-flop phenomena relies on the spin-spin 
coupling (fed by the relative spin-orbit precession in the 
unequal mass case), one can predict that not only BHs 


>, 

ro 

o 
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FIG. 22. Probability distribution of unequal mass flip-flops 
with respect to L-frame on top and with respect to J-frame 
below. These results are for initial random spin distributions 
evolved with 3.5PN form 100M of initial separation down to 
5 M (approximating the merger). 


can see their spins flipped, but material bodies with 
high spin, such as neutron stars, will as well, since they 
have been observed to pair with comparable masses. A 
simple evaluation of the flip-flop frequency and angles 
from Eqs. (19) and (30) applied to the double pulsar 
PSR J0737-3039 [67] with M = 2.587 M 0 , q = 0.935, 
M/r = 4.4 x 10~ 6 , leads to a flip-flop period of 1,230 
years, i.e., the spin direction of the less massive neu¬ 
tron star B would change by a degree every 3.4 years. 
However, the total change, according to Eq. (19) will be 
less than a degree if the highly spinning pulsar A has 
an intrinsic spin a-i < 0.25, i.e., in order to produce ob¬ 
servable effects nowadays its measured period should be 
below a millisecond instead of the actual 23 milliseconds. 
Notably, the spin precession frequency of pulsar B has 
actually been observed to be around 4.77° per year [651 . 
Unfortunately this precession of pulsar B meant that 
the highly beamed radio-emission is currently no longer 
pointing towards the Earth and we lost the detection 
of this one pulsar in 2008 69]. Yet the spin-orbit pre¬ 
cession of another pulsar has been measured. For the 
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PSR B1534+12 (a 37.9 milliseconds pulsar) a precession 
of 0.59° per year has been detected [TO]. This under¬ 
lines that other potentially observable neutron star bi¬ 
nary systems, at closer binary separations, may lead to 
larger observational flip-flop effects. 


3 q S 2 n S 
2 (1 + q)r 3 


1X -(! + <?) + 2 % 


where S 0 £ is a conserved quantity and 




i+¥) + vT <“) 
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Appendix A: Post-Newtonian Spin Evolution 
Equations in the L- and J- frames 

We discuss the spin evolutions in Eqs. ([7]) by decom¬ 
posing the spins along L and perpendicular to it, with 
unit vectors A and n, as shown in Fig. 
that this frame evolves as 


12 


Here, we have 


i 2 S 
L = -+V 


2 +i-ffn V Sq h ‘-’oL 


£ r 3 


A, 




(Al) 


\ ^ / 2 S e fffL 3 Tj Soft ^OL 

A - ^ 


L , 


where r? = q/(q + 1) 2 _, S eS = [1 + 3/(Aq )]§ 1 + [1 +3q/A]S 2 , 
and So = [1 + l/q\Si + [14- q]S 2 . Then, we obtain the 
spin evolutions from Eqs. Q as 


Sin — 


3 S, t S n ; 


TI+2A 

r ; 


— -fliSift + 

3S lL Sm 

qr 3 

• _ 3S m 

lL ~ 2 q(l + q)r 3 


3 qS lL S 2 n 


(*+§)• 


1 - 


2 S 0 l 


2 r 3 \ (1 + q)i 

<1 S ol 


1+ (1 + q)£, 

q(l + q)(2 + q)S 2 ~ x 


+2 ll + q + q^ ) S lX 


for Si and similar for S 2 exchanging labels 1 and 2 and 
q —>• 1/q. In the above expression Sn — SfL, Sin — Si'Tl , 
and S iX = Si ■ X; £ = \L\ = qM 3 l 2 r 1 l 2 /(I + q) 2 , and 


s l = s-l. 


1 . Flip-flop frequency 

The evolution of Si and S 2 along L becomes 
S lL = ^(Si x (S + 2S eS ))-L 

- 3 Soii -S lX (l + q S ° L 


r 3 1 + q 


1 


1 + q i 


S 2L = ^(S 2 x(S + 2S eB ))-L 


3 qSm S 2X fl + ^% ) • ( A4 ) 


r 3 1 + q 


1 + q l 


Then, the averaging over the orbital period (denoted by 
the brackets ()) gives 


<+> = Aa+«)Wi- W- 




3 I s ) I'Al (i Q S oL 


1L/ 2 r 3 ^ t y (1 + q ) 2 l 

x(5i x S 2 ) ■ J , 

^2b) = “2^3 ( 1 + “) — ~ (l + g )2 £ 

x(Si x S 2 ) ■ J. (A5) 

This shows the conservation of Sq ■ L, as we have (5j f ) + 

<l{S 2 i) = o. 

On the other hand, the spin evolution projected along 
J can be derived from Eqs. ([ 7 ]) as 

Su = = + (Si xS 2 )-L 


VI 


2 q 


+-^—es ofl s lX + ~^—S 0 n(Si xS 2 )-n) , 

1 + q 1 + q 


S 2 j - 




+ - 3 ^eS 0fl S 2 - x + ^-Son(S 2 X Si) ■ fl) , 

1 + 9 1 + q J 

(A6) 

(A2) and then the averaging over the orbital period gives 

(4j) = - +«) ( l ~ TTT^^) (Si xS 2 )-J. 


(1 + q ) 2 e 


(A7) 
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This shows the conservation of the quantity Sjy-J, where, 
Sn = qSi + S 2 , as we have q(S 1 f) + (S 2 j) = 0. 

Note that the factor £/(q\J\) (or qi/\J\) is different 
from (S x i) (or 0 2 i))- This difference can be expressed 
as 


(Su) ~ 



1 


\J\ 

1 

1 


0i • S - Si ■ t) 

0i -Si-Si- (j- &)) 


Si-S 2 


(A8) 


Therefore, the time evolution of the angle between Si and 
S 2 creates the different g-behavior into the two projection 
frames. This will affect the flip-flop angle, but not its 
frequency as we see below. 

To obtain the flip-flop frequency (measured along L or 
J will be the same) we first note that 


+ ^ s ^( 1 + YT ~ q f ¥) ’ 

0 oL )=O. (A9) 

Therefore, the time dependent part of the orbital- 
averaged first derivative of the spins, (S 1 i), (S 2 i), 0ij) 
and (S 2 j) is only (Si x S 2 ) ■ J, on which we focus in the 
following. 

The first derivative of (Si x S 2 )- J with respect to time 
can be written as 


>'■ * *' 


— 0 (1 _ 1 _ 

2 r 3 |j| V (1 + g) 2 t J 


(1 ~q 2 ) 


((L x Si) x S 2 ) ■ L 


L q 

- (s 0 X 01 x &)) • L 


(A10) 




1A 


01 x S 2 ) ■ L 


1 + 


9 Sql\ 

1 + 9 i J 


Here, we are always using the averaged evolution equa¬ 
tions.Taking one more time derivative of the above equa¬ 
tion, we have 


d0 

dt 2 


[0i x S 2 ) ■ J} 


_ 9 (1 — g) 2 (l + g) 2 l 2 + 9(l-q)(l + q)£S lL 
4 q 2 r 6 qr 6 

+ 9 (1-9 ? Sil s 2 l 9 (l-g) (5 + 3q)S 2t 2 

2 qr 6 4 r 6 

x0i x S 2 ) ■ J , 


9(1 - g)(l + q)iS 2L 9 (3 + 5 q) (1 - q) S lL 2 
qr 6 4 q 2 r 6 

9 (1 + 9) 2 (<Si 2 + 2 q Si ■ S 2 + q 2 S 2 2 ^ 

4 q 2 r® 

(All) 


where (-1—) means higher PN order terms. This equation has the form of an harmonic oscillator (with slowly varying 
frequency) and reaches the two extreme positions when Si, S 2 , and J (or equivalently at those points L) are all 
coplanar (and hence the triple vector product vanishes). 

Therefore, we can extract a frequency as 

_ 9 (1 — q ) 2 M 3 / | m | Q (1 -g)(Si L -S 2 L )M^ 9 (l-g) (3+ 5 g) S lL 2 
4 (1 + q ) 2 r 5 y r ) (1 + gjr 11 / 2 4 q 2 r 6 

9 (l-qfS lL S 2L 9 (1 — g) (5 + 3 g) S 2 0 9 Sq 

+ 2 g r 6 + 4 r® + 4 + "' ’ (A12j 


where we have used t = (g/(l + g) 2 )Af 3//2 r 1 / 2 (l + 2M/r) 
from Eq. (4.7) of Ref. [5T], and (+ • ••) stands for 
0(l/r 13 / 2 ) which is higher PN order corrections. Here, 
flff is constant at fixed r for g = 1. In the unequal 
mass cases, this is not constant in time, but we may 
treat it as a constant in the flip-flop time scale because 

oiSl ff /n ff = 0(l/rV*). 


2 . Maximum flip-flop angle 


The maximum flip-flop angle is calculated as follows. 
First, we obtain from Eq. (All I as 


0i x S 2 ) ■ J = Asin(Qfft ), 


(A13) 


where we have chosen the initial value by 5^(0) = — Si 
(cx,r (0) = — aq where a* = Si/rri 2 ) to obtain the maxi¬ 


mum flip-flop angle for the smaller hole. Using Eq. (A101, 
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the amplitude A is determined from 


x(S 2 2 -S 2L ( Of), 


(A 14) 


Sif , A= ±J-( 1 - 

ff 2r 3 |jj \ (1 + g) 2 £ 


q S ° L '(i + q )Si 


at t = 0 as 


A = 


(a 2 2 — a 2 i(0) 2S j a\ q 2 M 9//2 (a 2 2 — a 2 f( 0) 2 ^ (— di g 3 + 2 di g 2 + 2 d 2 p(0) q — d 2 p(0)) a\ qM 5 

(1 + qf (1 - q) \fr 

l (d 2 2 — a 2 i (0) 2 ^ aiM 11 / 2 


(i + qf (i ~ q ) 2 r 


\2 2 


(1 + q) (1 - g) 3 r 3 / 2 


(2di g 1 — 6«i g' — 6 a± q 4 a 2 f( 0) + 6di”g 4 + 20di q 3 a 2 f( 0) + 8d 2 £(0) g 


— 2 a 2 ~q — 6«i g 2 d 9 p(0) — 8ga 2 ^(0)“ + 2 a 2 q — d 2 + 3a 2 ^(0) ). 


(A15) 

where d 9 p(0) is derived later. Substituting this solution into Eqs. (A5) and integrating it with respect to time, we 
obtain 


T l 


t (t) _ («2 2 -« 2 i(0) 2 ) M 

r(l~qf 


ai 


(a 2 2 - a 2l (0 ) 2 ) (ot\ q + a 2i (0)) M 3 ' 2 

2 (! ~ cos (flfft)) + 2 A--A—^- (1 _ cos (flfft)) 


(1 — q) r 3 / 2 


[a 2 2 - a 2 f( 0) 2 ^ ^3di 2 g 2 + 10 a\ qa 2 f( 0) — a 2 2 + 4d 2 p(0) 2 ^ M 2 


r 2 (1 — q) 


(1 - cos (fi ff t)) 


(A16) 


with a 1 f(0) = — u\. Therefore, we have the maximized with the parameter d 2 p(0) is obtained as 
flip-flop angle at cos(f If ft) = —1, and d 2 p(0) is fixed to 
maximize the angle as 

a 2 ^(0) _ ol 2 \[M ' 3qa 2 a x M (A17) 


d\ 


+ 


a 2 (1 - q)^/r (1 - q) 2 r 

Then, the solution becomes 


a 


1 L 


tit) 


di 


= -1 + 


a 2 M 2qa 2 2 a\M 3 / 2 


+ 


(1 — q) 2 r (1 — g) 3 r 3 / 2 


x (1 — cos(U /ft)). (A18) 

As the result, the maximum flip-flop angle becomes 


a i L 0 max) 


oil 


+ 1 = 2 


a 2 2 M 2qa 2 2 a.\M 3 l 2 


(1 — q) 2 r (1 — g) 3 r 3 / 2 


(A19) 


Here, we need a special treatment for the q = 1 case 
because of the lack of the leading term in Eq. (A12). 


(a 2 2 - a 2t { 0 ) 2 ) 

ai 2 — 2 «i d 9 ^(0) + d 2 2 
(a 2£ (0)-ai) 2 
ai 2 — 2 a\ a 2 i(0) + d 2 2 


COS (flfft) 


(A20) 


The maximum flip-flop angle is found for a 2 ^(0) = di 
under an assumption oq < a 2 , and the solution becomes 




Oil 


= - cos (fi ff t). 


(A21) 


Therefore, the maximum flip-flop angle is tt. 

On the other hand, in the case of S 2 f(0) = S 2 
( a 2 f,(0) = a 2 )i w e can obtain the maximum flip-flop an¬ 
gle for the larger hole. The amplitude A is derived from 
Eq. dAl0b, 


riff a = 


3 £ 


1 - 


q 


Furthermore, since f Iff has only the leading order for q = 

1, we focus on the leading order calculation. Using the 

same initial configuration, a 1 ^(0) = — oq, the solution at t = 0 as 


2 qr 3 \J\ 
x^-S^O) 2 ). 




(A22) 


A = 


(di 2 - a 1 £(0) 2 ^ d 2 q 3 M 9 / 2 ^d! 2 - d 1 p(0) 2 ^ (d 1 | j (0) q 3 - 2 a 1 f(0) q 2 + 2 a 2 q - a 2 ) a 2 q 2 M 5 

(1 + qf (1 - q) i/r (1 + qf (i - qf r 

1 (ai 2 ~a 1 f j {£)) 2 \qa 2 M 11 / 2 

- - 4/1 ^ a3 ^ /9 -(3 a l£ (0) 2 q 6 - di 2 g 6 - 8d lZ/ (0) 2 g 5 + 2 di 2 g 5 + 6d-^O) g 4 d 2 + 8d lZ (0)'g 4 


) 2 g 6 _d 1 2g6_ 8 d li{0) 2g5 +2 ~- 2 ~5 , « ^ „4„._ , -/n+„4 

—2di 2 g 4 — 20d 1 ^(0) q J a 2 + 6d 1 | j (0) q 2 a 2 + 6a 2 2 q 2 — 6d 2 2 g + 2d 2 2 ), 


2 (1 + g) 4 (1 — g) 3 r 3 / 2 

2„4 on /nl «3 


where d 1 p(0) is derived later. From a similar analysis in 
the above, we have 


x (l-cos (rifft)), 


(A23) 

(A24) 


« 2 i(i) 


= 1 - 


d 2 


q 2 a\ 2 M 2g 2 di 2 d 2 M 3 / 2 


_(1 — g) 2 r (1 — g) 3 r 3 / 2 












































(A25) 


(A29) 


ot l i( 0) qai\fM 3qa\a 2 M 

ai (1 - q)y/r (1 - q) 1 2 r 


ai 2 (1 -q)y/r 

1 (q 2 + 4 q + l)aia 2 M 

2 (1 — q) 2 r 


Therefore, the maximum flip-flop angle is 


1 - 


Oi 2 L (^min) 


01 2 


= 2 


q 2 ot\ 2 M 
(1 — q) 2 r 


2 q 2 ai 2 a 2 M 3 4 5 6 ! 2 
(1 — g) 3 r 3 / 2 

(A26) 


In the J-frame, we find the maximum flip-flop of the 
larger hole with a 2 i (0) = a 2 and a parameter 0 ^( 0 ). 
The solution is 


a 2 


= 1 + 


q 3 ai 2 M q 2 (1 — 3q) ai 2 a 2 M 3 / 2 


(1-9) r 


(1 — q) r 3 / 2 


rr, 1 q ( q + X ) “1 M 

X COS- - ----2- 

2 (1 - g) r 

q [q 3 + 2 q — l) ai 2 a 2 M 3 / 2 
(1 — q ) 3 r 3 / 2 

and the maximum flip-flop angle becomes 


*2 jftmax) 
£*2 


l 2J 


1 (t min) 


r r»3 


= 2 


q d ai~M 
(1 — q) 2 r 
q 2 (l — 3q)ai 2 a- 2 M 3 1 2 


a 2 


(A27) 


(A28) 


(1 — q) 3 r 3 / 2 

To obtain this flip-flop angle, the initial configuration is 


a 2L(0) _ | 

ol 2 

a 2 j(0) ^ 1 q 2 a\ 2 M qa\ 2 ot 2 M 3 / 2 

a 2 2 r r 3 / 2 

« 1L (0) _ _1 (3 - g) gag \[M 
ai 2 (l-g),/r 

( 1 (g 2 — 8 q + l)oua 2 M 
+ 2 (1 -q) 2 r ’ 


For the smaller hole case, when we choose the initial 
value by £.^(0) = — S\ and find a 2 ^(0) to maximize the 
flip-flop angle that is the same ansatz as the L-frame 
analysis, the solution is derived as 


T j 


fit) 


Oil 


= -l+ - 


a 2 2 M (3 — q)a 2 2 otiM 3 / 2 \ 

q(l-qfr (1 - g) 3 r 3 / 2 J 


X COS (flfft) 


1 (q 2 + l) a 2 2 M 


2 q 2 (1 — q ) 2 r 
[q 3 — 2 q — l) a 2 2 a\ M 3 / 2 
q (1 — q ) 3 r 3 ! 2 

and the maximum flip-flop angle becomes 


Qyj(^max) Qyj(^mm) 


Oil 


+ 


= 2 

Ol\ 

(3 - q)a 2 2 a 1 


a - 2 2 AI 


72 


(1 — q) 3 r 3 / 2 


(A30) 


(A31) 


The above solution is obtained for the initial configura¬ 
tion, 


a lL(°) 

Oil 

Oil 

a 2i(°) 

Oi 2 


ol 2 j(Q) 

Ol 2 


- 1 , 

1 a 2 2 M a 2 2 ai M 3 ! 2 

— 1-1--- 1 —-—- - 

2 q 2 r qr 3 / 2 

1(1 — 3 q) a 2 \[M 
2 q(\-q)y/r 

1 ( q 2 — 8 q + l)a 2 aa M 

2 (1 — q) 2 r 

1 (1 + q) a 2 y/M 

2 q(\ — q) \Jr 

1 ( q 2 +4 q + l)a 2 on M 

”*"2 (1 — q) 2 r 


(A32) 
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